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FAST DOMINO TILEABILITY 


IGORPAK*, ADAM SHEFFERt, AND MARTIN TASSY* 


Abstract. Domino tileability is a classical problem in Discrete Geometry, famously solved 
by Thurston for simply connected regions in nearly linear time in the area. In this paper, we 
improve upon Thurston’s height function approach to a nearly linear time in the perimeter. 


1. INTRODUCTION 


Given a region R and a set of tiles T, decide whether R is tileable with copies of the tiles in T. 
This is a classical tileability problem, occupying a central stage in Discrete and Computational 
Geometry. For general sets of tiles, this is a foundational problem in Gomputability [Berl IBoas| 
and Gomputational Gomplexity |GJ[ILev[lVl] . For domino tiles, the problem is a special case of 
the Perfect Matching problem. It can be solved in polynomial time; even the counting problem 
can be solved at the cost of matrix multiplication (see e.g. |LPllK2] b 

In 1990, Thurston pioneered a new approach to the subject based on the study of height 
functions, which can be viewed as integer maps on the regions |Thu) . Thurston outlined a 
domino tileability algorithm which later has been carefully analyzed (see |5.4| ) and significantly 
extended to many other tileability problems (see (5.31. 


Theorem 1.1 (Thurston, 1990). Let R be a simply connected region in the plane 1?, and 
let n = \R\ be the area of R. There exists an algorithm that decides tileability of R in time 
O(nlogn). 


This result is worth comparing to the classical Hopcroft-Karp algorithm which has Oirt'/'^) 
time for testing whether a bipartite graph with n vertices and bounded degree has a perfect 
matching. This general bound has been significantly improved in recent years (see (5.1). 

Note that for polynomial time problems the cost of the algorithm depends heavily on how 
the input is presented. In case of graphs, the input is a list of vertices and edges, of size 0(n). 
On the other hand, the region R is traditionally presented as a list of squares, thus of size 
Ø(nlogn). 


The main idea behind this work is that plane regions R can be presented by the set of 
boundary squares. The input then has size Ø(plogp), where p = |9i?| is the perimeter of R. 
More economically, for simply connected R the input can be presented as a sequence of directed 
boundary edges Right, Up, Left and Down, starting at the origin. Of course, in this case 
the boundary set of squares can be computed in time Ø(plogp). Either way, one can ask if in 
this presentation one can improve upon Thurston’s algorithm. Here is our main result: 


Theorem 1.2. Let R be a simply connected region in the plane I?, and let p = |9ii| be the 
perimeter of R. There exists an algorithm that decides tileability of R in time O(plogp). 
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The result gives the first improvement over Thurston’s algorithm in 25 years, and is nearly 
optimal in this presentation. Clearly, the perimeter p = VL{^Jn) can be significantly smaller 
than n, so Theorem 1.2 improves upon Thurston’s for all regions with p = o(n). 


Let us note that Thurston’s algorithm not only decides domino tileability, but also constructs 
a domino tiling when the region is tileable. Specifically, Thurston shows that for every tileable 
region i? there is a unique maximum tiling To , corresponding to the maximum height function 
of the region R. He then inductively computes To in time O(nlogn). Clearly, it would be 
unhelpful to match this result, since listing all dominoes requires H(nlogu) time. However, we 
can do this in the following Oracle model. 

A tiling T of a region R is said to be site-computable with a query cost t if after preprocessing, 
for every square x G R, we can compute the adjacent square ?/ of a domino in T in time t. See §5.2| 
for the reasoning behind this model. 


Theorem 1.3. Let R be a simply conneeted region in the plane 1?, and let p = |i9i?| be the 
perimeter of R. If R is tileable with dominoes, the maximum tiling To is site-computable with 
a preprocessing time 0{p\ogp), and a query time O(log^p). 


The rest of the paper is structured as follows. In the next lengthy Section we present a 
criterion for domino tileability in terms of the height function on the boundary. The algorithm is 
presented in the following Sectionj^ We then show that similar results also hold for a triangular 
lattice (Section 1^. We conclude with final remarks and open problems in Section]^ 


2. Tileability theorems 

In this section we present necessary and sufficient conditions for the domino tileability of simply 
conneeted regions in Z^. In the following section, we use these conditions to construet an 
algorithm for cheeking the tileability of a simply conneeted region. For the rest of this section 
we assume that I? is endowed with a chessboard coloring, such that the square with corners 
(0,0) and (1,1) is white. 

Let i? be a simply conneeted region of Z^ such that the origin is on the boundary of R and 
Tfi is a tiling of R. A classic result (for example, see Fournier |Fouj 1 states that there exists a 
height function h : R —>■ Z that corresponds to Tr, and is defined as follows: 

• h{0, 0) = 0, and 

• for an edge {x, y) in Tr, such that when Crossing from x to y there is a white square to 
our left, we have h{x) — h{y) = 1. 

We denote by dR the points of Z^ that are on the boundary of R. To obtain a height function 
h : dR — >■ Z, we start from (0,0) and travel along the boundary. We begin by setting h{Q, 0) = 0 
and then place a height value to each point that we visit, according to the second condition 
of the definition. That is, when we cross an edge from point x to point y, if there is a white 
square to our left we set h{y) = h{x) — 1, and otherwise h{y) = h{x) + 1. If i? is not tileable, 
the final edge in our trip might not satisfy the condition. When this happens we say that dR 
has no valid height function. 

The following well known lemma characterizes the funetions h : i? —>■ Z that are height 
funetions of some tiling of R. 

Lemma 2.1 (see lEu])- Let R be a simply conneeted region of Z^ with the origin on its 
boundary, and let Jtfn be the set of height funetions that correspond to tilings of R. Then a 
funetion h : i? —)■ Z is in JtfR if and only if the following conditions hold. 

(i) For every x,y G R such that x = {ai,bi), y = ( 02 , 62 ), and oi — 02 = 61 — 62 = 0 
mod 2, we have h(x) — h{y) = 0 mod 4. 
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(ii) For every edge {x,y) G dR such that when Crossing from x to y there is a white square 
to our left, we have h{x) — h[y) = 1. 

(iii) For every edge {x,y) G R, we have \h{x) — h{y)\ < 3. 


One can also consider a tiling of the entire plane Lemma [2.l| remains valid for such 
tilings, with Condition (ii) becoming redundant. In this case, the set of height functions 
that are zero at a given point x = (a, b) has a maximum element a(x, •) (i.e., for every p G Z^, 
no height function h G satisfies h{p) > a{x,p)) that is defined as follows. For a point 
1 / G Z^, we write y = x + {i, j) and set S{i,j) = i — j mod 2. If a — & = 0 mod 2, then 


a{x,y) 

If a — 5 = 1 mod 2, then 

a{x,y) 


2||y-x||oo +S{iJ), 

if i > j, 

2||y-x||oo -S{i,j), 

if i < j. 

2||y-a:||oo -S{i,j), 

if * > j, 

2||y-x||oo +S{iJ), 

if * < j- 


It is easy to see that a{x, •) satisfies conditions (i) and (iii) of Lemma 2.1 and is thus a height 
function that corresponds to a tiling of 1?. Examples of maximal tilings that correspond to 
a{x, •) are depicted in Figurej^ (left and center). 




Figure 1. Left and center: Maximal tilings that are centered at x and corre¬ 
spond to the height function a{x, •). Right: Geodesic paths between x and y. 


We say that a sequence of points {xi ,..., a;„) is a geodesic path if 

• For every i < n, we have ||a:i+i - æiHoo = \\xi - xi||oo+l- 

• For every i < n, the points Xi and Xj+i are corners of a common 1x1 square in Z^. 


The right part of Figure [^depicts several geodesic paths between the same pair of points. 

One key observation in our proofs is that a{x, •) is strictly increasing on any geodesic path 
( æi , ..., Xn)- By the conditions of Lemma 2.1 for any 1 < i < n and height function h that is 
defined on Xi and Xi+i, there are exactly two possible values for h{xi+i) — h{xi). One of these 
two values is negative, and the other is ot{xi,Xi+i) — a{xi,Xi). That is, for any height function 
h we have h(xi+i) — h(xi) < a{xi,Xi+i) — a{xi,Xi). 

Another useful observation is that for every geodesic path {xi ,..., Xn) and 1 < i < n, we have 
that (xi,..., Xi) and {xi ,..., x„) are also geodesic paths. By combining this with the previous 
observation, we notice that a(x, •) is additive on geodesic paths. That is, if (xi,...,x„) is a 
geodesic path and 1 < i < n, then a{xi,Xi) + a{xi,Xn) = a{xi,Xn)- 

For a simply connected region R and two points x,y G R, we write x y if there exists a 
geodesic path between x and y that is fully in R. The following lemma gives a necessary and 
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sufScient condition for the tileability of a simply connected region R. This condition depends 
only on the height differences between pairs of points on the boundary dR. 

Lemma 2.2. Let R be a simply conneeted region of l? that contains the origin and let h : 
dR Z be a valid height function of dR. The region R is tileahle if and only if for every pair 
x,y G dR that satisfies x y, we have 

h{y) - h{x) < a{x,y). 


Proof. We first prove that the condition is necessary. Assume that R is tileable and extend the 
domain of h to all of R according to a specific tiling of R. Let (xi,a;„) be a geodesic path 
between two points x,y G dR (that is, Xi = x and Xn = y). Recall that for every 1 < i < n we 
have h{xi+i) — h(xi) < ot{x.,Xi+i) — a{x,Xi), which implies 

n—1 n—1 

h{y) - h{x) = ^ (h(xi+i) - h{xi)) < ^ {a{x,x^+l) - a{x,Xi)) = a{x,y) - a{x,x). 

i=l i=l 

Since a{x, x) = 0, we get that h{y) — h{x) < a{x, y) which completes this part of the proof. 

We next prove that the condition of the lemma is sufficient. For that, we show that the 
function 


( 1 ) 


x^oH , xr^jiy 


[h{x) + a{x,y)\ 


satisfies the three conditions of Lemma 2.1 (with respect to R). For Condition (ii), it suffices 
to show that for every y G dR we have hmaxiv) = h{y) (since h satisfies this condition by 
definition). Consider such a point y G dR, then h{y) + a{y,y) = h(y) and the assumptions of 
the theorem implies that for all x G dR such that x y the inequality h{x) + a{x,y) > h{y) 
holds. Hence hmax(y) = h{y) on dR, and Condition (ii) is satisfied by /i^ax- 

For every x G dR, the function h + a{x, •) satisfies Condition (i) since both h and a{x, •) 
are height functions. Since hmaxiy) = h{y) for every y G dR, the function hmax satisfies 
Condition (i) on dR. This “forces” the various functions a{x, •) to be identical mod 4 on dR, 
and thus all over R. That is, for any y G dR the expression h{x) + a(x, y) mod 4 does not 
depend on the choice oi x G dR. This in turn implies that /imax satisfies Condition (i). 

It remains to prove that hmax satisfies Condition (iii); that is, to show that for every pair 
of adjacent points x,y G R, we have |/imax(a;) — hniax(2/)| < 3. If both x and y are in dR, this 
is immediate from Condition (ii). Thus, without loss of generality, we assume that x is in the 
interior of R. Let z € dR satisfy h^s,x{x) = h{z)+a{z, x) and let P = (æi, ..., x„) be a geodesic 
path from z to x that is contained in R. Let Xi be the last vertex in P that is on dR, and 
notice that P' = (xi,..., Xn) is a geodesic path from Xi to x. Since hm&xix) has the maximum 
increase rate that any height function may have, we obtain that h^a,x{x) = h{xi) + a(xi,x). 
Since P' does not contain edges of dR, there must also exist a geodesic path from Xi to y. 
Since x and y are neighbors, this implies h^a_^{y) < h{xi) + a{xi, y) < h{xi) + a{xi, x) + 3 and 
hmaxiv) — h^ax{x) < 3. A symmetric argument yields hj^axix) — hjnaxiy) < 3, and completes 
the proof of Condition (iii). □ 


For points x,y G dR, we denote by G(x, y) the set of points in that are in at least one 
geodesic path between x and y. Notice that G{x,y) is a rectangle with edges of slopes ±1, 
possibly with two opposite corners truncated; for example, see the right part of Figurej^ Let 
S* be a set that contains dR and any number of points from the interior of R. We write x y 
ii x,y G S and G(x, y) \ {x, y} is disjoint from S. The following theorem is a refinement of 
Lemma [2.2[ which reduces the number of point pairs that determine whether a region is tileable. 

Theorem 2.3. Let R be a simply conneeted region oflf' that contains the origin, let h : dR -G Z 
be a valid height function, and let S G R be a subset that contains dR. The region R is tileable 
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Figure 2. If y is adjacent to the interior point x, then a geodesic path from 
x' to X implies a geodesic path of a similar length from x' to y. 


if and only if there exists a function g : S ^ Tj such that g = h on dR and for every pair 
x,y G S that satisfies x «s y, we have 

(2) - a{y,x) < g{y) - g{x) < a{x,y). 


Proof. We first prove that the condition is necessary. Assume that R has a tiling T and let 
g : i? —>■ Z be corresponding height function. By definition, g = h on dR and g satisfies ([^ for 
every x,y G dR with x y. For every pair x,y G dS for which x «s y there exists a geodesic 
path (xi,..., Xn) in R with x = xi and y = As in the proof of Lemma 2.2 we have 

n—1 n—1 

9{y) - gix) = ^ (gixi+i) - g{xi)) < ^ {a{x, x^+i) - a(x, Xi)) = a{x, y). 


i=l 


i=l 


A symmetric argument implies g{x) — g{y) < a(y,x), which completes the proof of this part. 

To prove that the condition of the theorem is sufficient, we show that it implies the condition 
of Lemma 2.2 That is, if a function g satishes Q for every pair x,y G dS with x «s y, then 
the same condition is also satisfied for every pair x,y G dR with x y. 

Consider a pair x,y G S such that x '^r y. We prove that ([^ holds for a ;,y by induction 
on \\x — y\\oo- Since dR C S, this would complete the proof of the theorem. For the induction 
basis, consider the case where ||a; — y\\oo = 1- In this case we have x «s U, so ([^ is satisfied for 
X, y by the definition of g. 

For the induction step, consider the case where ||a; —y||oo = fc > 1. In this case, either x «5 y 
and (§ is satisfied by the definition of g, or there exists a geodesic path P = [xi, ..., x„) between 
X and y that is in R and intersects S\{x, y}. In the latter case, let Xi be a vertex of P that is in 
S \ {x, y}. Then (xi,..., Xi) is a geodesic path between x and Xi and {xi, ..., x„) is a geodesic 
path between Xi and x„. By the induction hypothesis, we have 

x) < g{xi) - g{x) < a{x, Xi), 
a{y, Xi) < g{y) - g{xf) < a{x^, y). 


By combining these two inequalities we get g{y) — g{x) < a{x,Xi) + a{xi,y). Since P is 
obtained by combining a geodesic path from x to Xi together with a geodesic path from Xi to y, 
we have a{x, Xi) + a{xi,y) = a{x, y), which in turn implies g{y) — g{x) < a{x, y). A symmetric 
argument implies g{x) — g{y) < a{y,x), which completes the induction step and the proof of 
the theorem. □ 


3. Algorithm for tileability 


3.1. Outline. In this section we prove theorems 1.2 and |I.3[ First, we present an algorithm for 
checking whether a simply connected region R is tileable. The algorithm is based on partitioning 
R into interior-disjoint squares of various sizes. These squares have their vertices in Z^, but are 
“rotated by 45° ” in the sense that the slopes of their edges are ±1. To cover R with such squares. 
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along the boundary of R we use right-angled triangles with two edges of length 1, instead of 
squares; for example, see Figure We consider the set S that consists of dR together with 
the vertices of the rotated squares. By Theorem |2.3[ to check whether R is tileable it suffices 
to compare between pairs x,y G S that satisfy x ~s D- We will prove that each point of S 
participates in at most eight such pairs, which would in turn imply that total number of pairs 
that satisfy a: «5 y is at most linear in the perimeter of R. 



Figure 3. A subdivision of the area that is bounded by the solid edges into 
rotated squares and right-angled triangles with two edges of length 1 . 


3.2. Partitioning the region. We begin with the following technical result. 

Theorem 3.1. Let R he a simply connected region with |5i?| =p. Then there exists a subdivi¬ 
sion S of R into 0{p) interior-disjoint rotated squares and right-angled triangles with two edges 
of length 1. Such a subdivision can be found in 0{p\ogp) time. 

Proof. We prove the theorem by presenting an algorithm that receives a simply connected region 
R with perimeter p and constructs a subdivision S' of i? in time 0{plogp). The subdivision S 
consists of 0{p) rotated squares and right-angled triangles with two edges of length 1. All of 
the squares in this proof are rotated by 45°. 

We begin the algorithm by initializing several variables. Set Qin = Qout = 0 where Qin 
(respectively, Qout) is a set in which we place squares that are fully on the inside of R (resp., 
fully on the outside of R). Let n be the smallest power of 2 that is larger or equal to 2p, consider 
an \/2n x square that fully contains R, and let Sq be a set that contains only this square. 
We repeat the following process for t = log 2 n iterations: 

• At the beginning of iteration i we consider the set which contains interior-disjoint 

squares of size x that were obtained in the preceding iteration. We 

partition each of these squares into four interior-disjoint squares of size •\/2 ?t-/2* x •\/2 R-/2b 
Denote the set of these 4|5'i_i| squares as R. For each square of we record the four squares 
of Qin U Qout U Ti that share a boundary with it (some of these might not be squares but the 
area outside of the square of S'o). This can be done in constant time by using the information 
that was stored in the previous iteration for the squares of Ti_i. 

• We travel across dR. Every time that we enter a square s € Ti, we mark the point from 
which we entered s, mark what side of this boundary point of s is the interior of R, and insert 
s into Si. 

• For every square s G Ti \ Si, we check whether s is in the interior or in the exterior of R (see 
below for the full details of this process). If s is in the interior, we add s to Qin. Otherwise, 
we add it to Qout- 

After log 2 n iterations, we have a set Qin of interior-disjoint squares that are fully contained 
in R and a set Siog^ n of x squares whose interior is intersected by dR. We split every 
square s G Siog^n into four right-angled triangles and insert into S the triangles that are in 
the interior of R (out of the four triangles, between one and three are in the interior). After 
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also inserting Qin into S, the set 5 is a subdivision of R into interior-disjoint squares and 
right-angled triangles. 

We now explain how, at the end of the i-th iteration, we go over each square s G Ti\Si and 
check whether s is in the interior or in the exterior of R. We go over the squares of R \ Si in 
an arbitrary order. When considering a square s G Ti \ Si, we already know which squares of 
Qin U Qout U Ti share a boundary with s. Notice that there exists a unique square in each side, 
and that these four squares may be of different sizes. 

• If One of the four surrounding squares is in Qin, we add s to Qin- 

• Otherwise, if one of the four surrounding squares is in Qoutj we add s to Qout- 

• Otherwise, if one of the four neighboring squares s' is in Si, we travel along the boundary of 
s' until we get to an intersection with the border of R (we marked these intersection points 
when inserting s' to Si). For each such intersection point we previously marked which side is 
the interior of R, and we can use this information to determine whether s is on the outside 
or on the inside of R (and then place s accordingly in or in Qout)- 

• We remain with the case where the four neighboring squares are currently in R \ {Si U Qin U 
Qout)- In this case, we arbitrarily choose one of these four neighbours s' and add s to the 
“waiting list” of s' (see below for the purpose of this list). 

If several squares of Ti \ Si form a connected component, then either all of these squares are 
in the interior of R or all of these square are in the exterior of R. Thus, each time that we 
decide whether a square s G Ti \ Si goes into Qin or into Qout, we inspect the waiting list of s 
and place the squares that are in it in the same Q (we then have to check the waiting lists of 
each of these squares, and so on). 

The running time of the algorithm. Notice that 2p < n < 4p, so for any asymptotic bound 
that we derive with respect to n, we may replace n with p. To show that the running time of the 
algorithm is 0{p\ogp), we require following lemma. Recall that Si is the set of interior-disjoint 
squares of size x at step i whose interior is intersected by dR. 

Lemma 3.2. \Si\ < 9 • 2*“^. 

Proof. Partition the square of Sq into 4® interior-disjoint squares of size yfinjT' x ^J^njT', and 
denote the set of these squares as S[. Notice that Si C S{. Specifically, Si consists of the squares 
of S[ that are intersected by dR. We traverse dR starting from an arbitrary point vi. During 
this process we will mark fewer than 9 • 2®“^ squares, so that the marked squares fully contain 
the boundary of R. This would immediately imply |S'i| < 9 • 2®“^. 

We first mark a square of S{ that contains vi (there are at most four such squares), and the 
eight squares that surround it (i.e., share a vertex with it). We then continue to travel across 
the boundary of R until we get to a point V 2 that is not contained in any marked square. We 
mark a square of Sl that contains V 2 and the eight squares surrounding it (some of these squares 
are already marked, and remain so). We then continue to travel until we reach a point that 
is not in any marked square. We continue in the same manner until we return to Vi. 

Notice that each time that we get to a point Vi that is in no marked square, we mark fewer 
than nine unmarked squares. After marking these squares, we travel at least 2n/2® steps along 
dR before we reach Vi+i. Since dR is of length p < n, the total number of marked squares is 
smaller than 9 • 2®“^ □ 


The algorithm consists of t = log 2 n iterations. Let us show that each iteration has a running 
time of 0(n), which would complete the proof of the theorem. Consider the running time of 
the i-th iteration. By lemma 3.2 we start this iteration with a set Si-i of 0(2®) squares, and 


partition it into a set R of 0(2®) squares. For each new square we also record the four squares of 
Qin U Qout U Ti that are its direct neighbors. Since handling each square of R requires constant 
time, this step takes 0(2®) = 0(n) time. 
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We then travel the boundary of i?, and every time that we cross to a different square of Ti 
we perform a constant number of operations (marking the entry point and possibly inserting 
the square into Si). By considering the origin to be the bottom left corner of the square of 
S'o, we can easily decide when we enter a new square of Ti. This occurs exactly when the x or 
y coordinate of our current position becomes 0 mod (njT). Thus, the entire traversal of the 
boundary of R takes 0{jx) time. 

The last part of the i-th iteration involves going over each square s € Ti \ Si and checking 
whether s is fully in the interior or in the exterior of R. This check is based on the the four 
squares that surround s. The only case that takes more than a constant time occurs when none 
of these squares is in and Qoutj while at least one is in Si. In this case we travel along 
the boundary of such a neighboring square. The peri mete r of such a square is 4n/2\ so each 
instance of this case takes OinjT) times. By lemma 3.2 |S'i| = 0(2*) and we consider each 


square of Si at most four times (at most once for each of its four direct neighbors). Thus, the 
combined time of all of these checks is 0{n). 

The only issue that we did not consider so far is the time required to handle the waiting 
lists. Since each square of Ti is in at most one such list, and since \Ti\ = 0(2*), the total time 
for handling the waiting lists is 0(2*) = 0(n). Finally, It is easy to see that the last step of the 
algorithm, of cutting the squares of St into triangles, requires 0(n) time. This completes the 
proof. 


•ti-l 


Bounding |S'|. By Lemma 3.2 in the i-th iteration the algorithm adds fewer than < 9- 
squares to Qin. Summing this quantity over the log 2 (n) iterations of the algorithm yields 
l^l = 0(n) = 0(p), as asserted. Some of these squares may be split into two triangles, but this 
does not affect the asymptotic size of S. □ 


Proof of Theorem \1.S\ We begin by running the algorithm of Theorem |3.1[ to obtain a subset 
S C Rof 0{p) interior-disjoint squares and triangles that cover R. We build a graph G = {V, E), 
where V consists of the points of dR and the vertices of the squares and triangles of S. An 
edge {x,y) S is in E if and only if x «s y. Notice that \V\ = 0{p). 


Bounded degrees. We now prove that every vertex of G is of degree at most eight. That is, 
that any point x G S satisfies x «s y for at most eight points y G S \ {æ}. 





'*( 


(b) 


Figure 4. (a) When x G dR, it is blocked by triangles. (b) The point x is the 
corner of three squares and one pair of triangles of S. 
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We first consider the case where x € dR, which forces x to be adjacent to triangles of S 
(and to at most one square); for example, see Figure|^a). We say that x forms a valid pair 
with y G V if X ~s U- In the current case, x forms a valid pair with vertices that are are in 
a common triangle with it. Moreover, x forms a valid pair with a point y G V that does not 
share a triangle with it if and only if the straight-line segment between x and y is fully in the 
interior of R, does not contain any other points of F, and has a slope ±1. Thus, in this case x 
participates in at most seven valid pairs. 

Next, consider the case where x ^ dR and is at the corner of four squares and/or pairs of 
triangles of S\ for example, see Figure|^b). As before, x creates a valid pair with each vertex of 
V that shares a triangle with x. The maximum degree of eight is obtained when x is surrounded 
by four pairs of triangles. If a; is a vertex of a square s G S, denote the two edges of s that are 
adjacent to x as ei and 62 - Notice that x creates a valid pair with the point of V that is closest 
to it along ei and with the point of V that is closest to it along V 2 ■ For example, in the case of 
the square to the right of x in Figurej^b), these vertices are yi and y 2 - 

Beyond the valid pairs that are described in the previous paragraph, x cannot create a valid 
pair with any other point of V. For example, in Figure |^b) x cannot form a valid pair with 
any additional vertex 2 : that is to its right, since there must be a geodesic path between x and 2 
that contains either yi or 1/2 ■ 



Figure 5. The case where x is on the boundary of a square s G S without 
being a vertex of s. 


Finally, it is possible that x ^ dR and x is on the boundary of a square s G S without being 
a vertex of s; for example, see Figure In the figure, x does not form a valid pair with 22, 
23, and 24, due to geodesic paths that contain 24. While x «s y2, we have x 9^5 yi due to a 
geodesic path that contains 2 / 2 - Similarly, x 765 j /3 due to x'. By the way in which we perform 
our subdivision of R into squares, it is impossible to have a subdivision with x but without x' 
(that is, the square below x cannot exist without the square to the right of x. The latter square 
may be further subdivided). Similarly, x' forms a valid pair with 1/3 but not with y 2 , j/i, and 
the vertices to the right of 2/3. 

The examples in the previous paragraph illustrate a general principle: When x is on the 
boundary of a square s G S without being a vertex of s, out of the points y gV ior which the 
segment xy intersects the interior of s, at most one point creates a valid pair with x. Specifically, 
such a point y creates a valid pair with x if and only if the segment xy has slope ± 1 , does not 
contain any other point of S', and is fully in R. For any other such point z G V, there must 
be a geodesic path between x and 2 that passes either through y or through one of the two 
neighbors of x along the boundary of s. Thus, in this case x is of degree at most six (and this 
degree is obtained when x is adjacent to two pairs of right-angled triangles). 
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Computing E. By the above degree restriction, we have \E\ = 0{p). To build E, we go over 
each vertex oi x & V and look for the other points of V that form a valid pair with x. By 
considering the above cases, we notice that if x forms a valid pair with y G V, then either x 
and y are on a common triangle of S or the segment xy has a slope of ±1 and no other points 
of V on it. To handle the former case, we simply go over every triangle in S and add its three 
edges to E. 

For every line i of slope 1 that contains points of V, we keep an array of the points of V 
that are on £, sorted by their æ-coordinate. There are 0(p) lines with a total of 0(p) points 
on them. Thus, the arrays can be built in O(plogp) time. We then go over every array and 
add an edge between every two adjacent points on it, with the following exception. If we get 
to a point x € dR on we check whether t leaves R in x and if so do not add the edge that 
intersects the outside of R. We repeat the same process for lines with slope -1, which completes 
the construction of E. Notice that this construction takes 0{p\ogp) time. 

By Theorem 2.3 i? is tileable if and only if there exists a height function g ■. S ^ 1 that 
satisfies ([^ for every x,y £ S with {x,y) £ E. We now describe an algorithm for finding such 
a function g (or stating that such a function does not exist) in 0{plogp) time. Specifically, 
out of the set of functions that satisfy the above condition we find the maximum function hmax 
defined in Q. 


Computing g. We begin the algorithm by initializing several variables. Let h : dR -£ Z he 
a valid height function. The beginning of Section explains how to hnd such a function in 
0{p) time (if such a function does not exist, we stop the algorithm and announce that R is not 
tileable). Let A be an array with a cell for every point x £ S, such that eventually we would 
have A[x\ = g{x) (that is, A would describe g). We initially place “N” in each cell of Al, to 
state that g is currently undefined for the point that corresponds to the cell. Then, for every 
X £ dR we set A[x\ = h{x). 

Let H he a heap with an element for every point of S (a standard binary heap would suffice). 
For each x £ S with A[x\ = N, the key H[x\ of x is the maximum integer value for g{x) that 
does not violate ([^ with points that already have a value in A. At hrst, we insert every point 
of S' to iL with key E[[x\ = oo. Every time that we update a cell A[x\ (including during the 
above initialization of the points of dR in A), we remove x from El and update the keys of each 
y £ S that is adjacent to a: in G (that is, for which x «s y)] specifically, for every such y we set 
H[y\ = iLnin{H[y], A[x] + a{x,y)}. Notice that setting a value in a cell of A leads to updating 
at most eight elements of H. 

The main part of the algorithm consists of repeating the following process until the heap H 
is empty and no cell of A contains the undefined value N. Let x be the point with the smallest 
key k in H. We set A\x\ = fc, remove x from H, and update the keys of points that are adjacent 
to a; in G as described above. We then go over the vertices that are adjacent to a: in G and 
already have values in A. For each such vertex y, we check whether x and y satisfy (i). If not, 
then we stop the algorithm and announce that R is not tileable. 

If the above process ended since the heap H is empty and no cell of A contains the undefined 
value N, then we obtain a function g that satisfies ([^. In this case, we announce that R is 
tileable. 


Correctness. To prove that the algorithm is correct, it suffices to prove that the function g 
that the algorithm computes is indeed the maximum height function hmax from 0 (although 
defined only on the points of S). By Theorem 2.3 R is tileable if and only if such a function 
exists. 

We first claim that for each y £ S \ dR, there exists a point x £ S such that x «s y and 
g{y) — g{x) = a{x, y). Indeed, if no such x exists then the algorithm would have assigned a larger 
value to g{y) after removing y from H. Since a{x, y) > 0 when x ^ y, we obtain that g has no 
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local minimum (with respect to the edges of G) outside of dR. Specifically, a straightforward 
induction on the (edge) distance of x from dR shows that for every y in the interior of R there 
is a path (xi, ... ,Xn) in G such that y = xi, Xn S dR, and g{xi+i) — g(xi) = a{xi+i,Xi) for 
every 1 < i < n. 

We show that g = /imax on S by induction on the number of values that the algorithm 
already set in A. For the induction basis, the claim follows by definition for every point of dR. 
For the induction step, assume that in the i-th iteration of the algorithm the point y G S \ dR 
is chosen, since it has a minimum key in H . By the induction hypothesis, this key is 

H[y] = min{A[x]+a{x,y)) = min + a{x, y)), 

where the minima are over all x s.t. {x, y) G E and A[x\ ^ N. 

From above, there exists a path (cci,... ,x„) in G such that y = xi, € dR, and g[xi+i) — 
g(xi) = a(xi_|_i,Xi) for every 1 < i < n. That is, H[y] = h{xn) + a(x„,y). By inspecting the 
definition of hmax in 0 , we notice that to have /imax( 2 /) = h(x„) + a(x„, y) it remains to prove 
that there is no point z G dR such that h{xn) + a{xn,y) > h(z) +a{z,y). A priori, this can only 
happen if the value h{z) + a{z, y) was not yet discovered by the algorithm since some geodesic 
path between y and z contains a point z' G S that is still in the heap H. 

To see why the problematic scenario is impossible, notice that in the i-th iteration the key 
of every point x that is still in H is at least h(x„) + q;(x„, y); otherwise we would have removed 
X from H before removing y. While some of these keys may be decreased in following steps of 
the algorithm, no key will be decreased to a value that is smaller than /i(x„) + a(x„,j/) + 1 . 
Thus, when removing z' from H, we have H[z'] > h{xn) + a{xn,y). This in turn implies 
h{z) + a{z, y) > H{z') + a{z', y) > h{xn) + Q!(x„, y). That is, the problematic scenario cannot 
occur, and the correctness proof is complete. 

Running time of the algorithm. As already mentioned, obtaining S requires O(plogp) 
time, building G requires O(plogp) time, and h is computed in 0{p) time. Computing g (that 
is, computing A) requires 0{p) steps. Every step involves a constant number of operations. 
Most of these operations require a constant time, except for removing the minimum element 
from the heap and updating the keys of at most eight other elements. Each such operation 
takes O(logn) time, so computing g requires O(plogp) time. □ 


Proof of Theorem [H We begin by describing the preprocessing step. In this step, we first 
run the algorithm of Theorem 1.2 to obtain a subset S C R oi vertices of 0{p) interior-disjoint 
squares and triangles that cover R. We also obtain the values of the maximum height function 
for the points of S. As stated in Theorem 1.2 this can be done in 0{p\ogp) time. We then 
preprocess the subdivision of R for point location queries (see e.g. |EGS| 1. Specifically, after 
a preprocessing time of 0{p), for any point x G R we can find the square or triangle that 
contains x in O(logp) time. 

We now move to describe the query step, where we are given a query point x G R. First, 
we consider the case where x is not on the boundary of any square in the subdivision of R. By 
using the point location algorithm, we find the square that contains x in O(logp) time. We then 
partition this square into four subsquares of equal size, and add the vertices of these subsquares 
into S. There are at most five new vertices, and we compute the height value of each in O(logp) 
time, as in the proof of Theorem 1.2 (using the 0{p) arrays that we built in that proof). Out 
of the four subsquares, we find the one that contains x (again, assuming that x is not on the 
boundary of any of them) and subdivide it into four squares as before. We repeat this process 
until X is surrounded by eight vertices of S, as in Figurej^ We then have the height values of 
X and of the eight vertices that surround it. By the conditions of Lemma |2.1| this implies the 
behavior of the maximum tiling of R around x. 
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Figure 6. The edges of the square Q are solid, while the edges that are added 
in the subdivision steps of the query are dashed. 

If at some point during the above process x is on the boundary of a square, we continue 
the subdivision process with two other points — the one immediately above x and the one 
immediately below x. Notice that such a split from x into two other points can occur at most 
once, and that at the end of the process we still obtain the height values of x and the eight 
vertices that surround it. 

In summary, the query algorithm stops after O(logp) steps, each taking O(logp) time. Thus, 
handling a query requires O(log^p) time. □ 

Remark 3.3. We believe that with slight modihcations the running time of the query can be 
improved to O(logp), but chose not to pursue this direction at this point. 


4. Lozenge tilings 

The results of this paper can be extended to other lattices, and specifically to lozenge tilings in 
the triangular grid, which are dual to perfect matchings in a hexagonal grid. In this section we 
present a brief outline for how to extend our result to the case of lozenge tilings. Once again 
we follow Thurston, who defined the corresponding height function in the original paper |Thuj 
(see also [ChallRl] ). 



Figure 7. Coloring the triangular grid with three colors. 

Let T be the triangular grid with a fixed coloring of the vertices in black, red, and blue, such 
that every edge is adjacent to two vertices of different colors (see Figurej^. Let vi be a vector of 
length 1 in the positive direction of the x-axis, let V 2 be vi rotated counterclockwise by 120°, and 
let V 3 be vi rotated clockwise by 120°. Every point of T can be written as avi + bv 2 + cv^ (where 
a, 5, c are non-negative integers) in infinitely many ways. For example, the origin can be written 
as (a, o,a) for every non-negative a. However, after adding the condition min{a, &, c} = 0, for 
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every point p of T there is a unique way of writing p as avi + bv 2 + CV 3 . Using these unique 
values of a, b, c, we say that the coordinates of p are (a, &, c). 

Let i? be a simply connected region in T. Similarly to the characterization of height functions 
for domino tilings in Lemma |2.1[ a function : i? —>■ Z is the height function of a lozenge tiling 
if and only if: 

(i) For every two vertices x,y G R, we have h(x) — h(y) = 0 mod 3 if and only if x and y 
have the same color. 

(ii) For every edge {x, y) G dR with respective colors (black,red), (red,blue), or (blue,black), 
we have h(x) — h{y) = 1 . 

(iii) For every edge {x,y) G i?, we have \h{x) — h{y)\ < 2. 

Also as in the case of domino tilings, for any x G R there exists a maximal height function 
of the plane ot{x,-) with a(x,x) = 0. Specifically, a{x,y) is defined as the sum of the three 
coordinates of y when considering x as the origin. In the triangular grid, we say that a sequence 
of points {xi,..., Xn) is a geodesic path if 

• For every i < n, we have Ijaii+i — a;i||i = \\xi — cciUi+l. 

• For every i < n, the points Xi and Xi+i are corners of a common 1x1 triangle in T. 

Unlike the case of Z^, in the triangular grid geodesic paths do not necessarily minimize 

the number of edges. Moreover, when travelling along a geodesic path, we only move in the 
directions 111 , 112,^3 (since moving in one of the other three directions will result in a step of 
distance 2 ). Specifically, every geodesic path uses at most two of these three directions. As 
before, for x,y G T we denote by G{x, y) the union of the geodesic paths between x and y. It 
is not difhcult to verify that G{x,y) is always a parallelogram (possibly of width zero). Let S 
be a set that contains dR and any number of points from the interior of R. We write x «s y 
when x,y G S and G{x, y) \ {cc, y} is disjoint from S. The following theorem is the lozenge tiling 
analogue of Theorem |2.3[ 


Theorem 4.1. Let R be a simply connected region in T that contains the origin and let S G R 
be a set that contains dR. Then R is tileable if and only if there exists g S ^ 'L such that 
g = h on dR and for every pair x,y G dS with x pc s U we have 

(3) - a{y,x) < g{y) - g{x) < a{x,y). 


The proof of Theorem |4.1| follows verbatim the proof of Theorem |2.3| We omit the details. 
Similarly, we can use a variant of Theorem 3.1 to find a subdivision of R into 0{p) interior- 
disjoint triangles of various sizes. The main difference is that here we use equilateral triangles 
with side length 2* instead of squares. Indeed, every such triangle can be subdivided into four 
interior-disjoint triangles of side length 2®“^. We set S to consist of dR together with the 
vertices of the triangles of the subdivision. For any x G S, it is not difhcult to show that at 
most six points y G S for which x pc s U (the actual bound seems to be smaller than six, but 
this does not matter for our purpose). Finally, by revising the algorithm of Theorem 1.2 we 
obtain the following theorem. 


Theorem 4.2. Let R be a simply connected region in the triangular grid of the plane, and let 
p = |9i?| be the perimeter of R. Then there exists an algorithm that decides tileability of R in 
time O {p logp). 


5. Final remarks and open problems 

5.1. For general bipartite planar graphs, recent developments improve the Hopcroft-Karp 
bound |HKj to nearly linear time. First, the existence of a perfect matching is equivalent to the 
circulation problem, where all white vertices have supply 1 and all black vertices have demand 1. 
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It is known that the circulation problem in planar graphs can be solved within the same time 
bound as the shortest path problem with negative weights on a related planar graph |MNj . The 
latter can be solved in time 0(n log^ n/log logn), see [MW| . This almost matches Thurston’s 
original O(nlogn) bound in Theorem o 

5.2. Our Oracle model for the perfect matching is similar to other models of sparse graph pre- 
sentations, which are popular in the study of graph properties of massive graphs (see e.g. [Gdl]). 
The idea is to give a sublinear size presentation of a perfect matching, amenable to running 
further sublinear time algorithms; see e.g. [RSj for a primer on the subject. 

5.3. Thurston originally defined and studied height functions for the domino tilings as in this 
paper and for the lozenge tilings in a triangular lattice. Since then, many generalizations and 
variations have been discovered. These include other rectangles in the plane [KKl IKoml IR2| . 
other tiles in the triangular lattice [Rlj . rhombus tilings in higher dimension |LMNj . perfect 
matchings of more general graphs in the plane and other surfaces |Chal lltol ISTCR j , and even 
infinite domino tilings [RFR] . We refer to [Pik] for a (somewhat dated) survey of various 
tileability applications of height functions and tiling groups. 

On the complexity side, there are a number of N P-completeness results for the decision and 
counting problems for general regions with small tiles, see e.g. |BNRRl IMR) . and more recently 
for simply connected regions [PYlj . In case of domino tilings, there are also # P-completeness 
results for 3-dimensional regions |PY21 IV2) . 

The notion of height functions for domino and lozenge tilings has also made a remarkable 
impact in Probability and MCMC studies (see e.g. |LRS[ lK2] b 

5.4. The complexity of Thurston’s algorithm has been investigated to a remarkable degree 

in the Computational Geometry literature. These include generalizations to regions with 
hoies parallel computing |Fou| . and more general graphs [Cha]. 

The idea behind our tileability criterion, stated in Lemma |2.2[ was first given in the third 
author’s thesis [T^ . in the context of tromino tilings. The criterion is especially surprising 
given the fundamentally non-local property of the domino tileability, as elucidated by the aug- 
mentablity problem (see [Korni §11.3]). 

5.5. We believe that our approach can be further extended to a variety of tiling problems 
which admit height functions, such as tilings with bars (see [BNRRI IKKI ITas] l. In a different 
direction, the heart of the proof is the idea of scaling represented by the squares which are used 
heavily in Section[^ It would be nice to see this idea can be further developed. Finally, if the 
boundary dR is given by some kind of periodic conditions (cf. mi), one can perhaps further 
speed up the domino tileability testing. Unfortunately, at the moment, we do not know how to 
formalize this problem. 

Acknowledgements. We are very grateful to Scott Garrabrant and Yahav Nussbaum for in- 
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